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We derive the exact solution of a one-dimensional Markov functional model with log-normally dis- 
tributed interest rates and constant volatility in discrete time. The model is shown to have two 
distinct limiting states, corresponding to small and asymptotically large volatilities, respectively. 
These volatility regimes are separated by a phase transition at some critical value of the volatility, 
at which certain expectation values display non-analytical behavior as a function of volatility. We 
investigate the conditions under which this phase transition occurs, and show that it is related to 
the position of the zeros of an appropriately defined generating function in the complex plane, in 
analogy with the Lee- Yang theory of the phase transitions in condensed matter physics. 
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I. INTRODUCTION 

An important class of interest rate models, which in- 
cludes many of the models currently used in practice, is 
the class of Markov- functional models [iHS] . The advan- 
tage of these models is that the value of discount bonds 
can be expressed as a functional of a low-dimensional 
Markov process. The specification of this functional de- 
pendence allows one to model the distribution of the 
forward rates with a prescribed probability distribution. 
This implies that these models can be calibrated exactly 
to a set of market instruments, such as caplets or swap- 
tions. 

We consider a one-dimensional Markov-functional 
model with log-normally distributed forward Libors in 
discrete time. We show that the model can be solved ex- 
actly for a constant (time-independent) Libor volatility, 
and exact results can be found for the dependence of all 
discount bonds on the Markovian driver. 

The exact solution of the model is used to study its be- 
havior as a function of volatility. We show that the model 
has two distinct regimes, in the low and large volatility 
limits, respectively. These regimes are separated by a 
sharp transition, occurring at some critical value of the 
volatility. We investigate the nature of this transition, 
and discuss the conditions under which it occurs. This 
volatility transition is similar to a first order phase tran- 
sition in condensed matter physics, and is described by 
an analog of the Lee- Yang theory of phase transitions 



II. THE MODEL DEFINITION 

We consider a Markov functional model with discrete 
time evolution. The tenor structure is a finite set of dates 



^ to < ti < ■ ■ ■ < t„ 



(1) 



representing maturities equally spaced, e.g. by 3 or 6 
months apart. 



The fundamental dynamical quantities are the zero 
coupon bonds Pij = Pti,tj- They are functions of a 
one-dimensional Markov process a:(i), which will be as- 
sumed to be a simple Brownian motion with the usual 
properties E[a;(i)] = 0,E[a;^(t)] = t. The model is de- 
fined by the probability distribution of the forward Libor 
rates Liiti) = ^{P~l^^ - 1) for the (ti,ii+i) period, with 
Ti = ti+i ~ti . We will work throughout in the i„— forward 
measure, with numeraire the discount bond Pt,tn- Specif- 
ically, in this measure the Libor rates Li will be assumed 
to be log-normally distributed 



L,, 



1 2 

Li exp ( ipXi - -ip t, 



(2) 



For notational simplicity we denote the value of the 
Markov driver at time ti as Xi = x(ti). We assume that 
the Libor volatility '0 is a constant, although a more gen- 
eral formulation of the model is possible, wherein tp has 
term structure. For the purpose of illustrating the phe- 
nomenon considered here, it will be sufficient to consider 
a constant Libor volatility ip. 

We denoted in Eq. ^ with Li the convexity-adjusted 
Libors; they are the expectation values of the Libor rates 
in the measure considered, and their determination is 
part of the solution of the model to be discussed below. 

The continuous time limit of this model follows the 
short rate process 



dr{t) 



iljdx{t) + {— logf {t))dt 



(3) 



where we introduced the convexity-adjusted forward 
short rate f{t) as the continuous time analog of Li. The 
model Eq. ([3]) describes a simple log-normally distributed 
short rate model without mean reversion [y, |7| . As dis- 
cussed in Ref. p| , mean reversion can be introduced by an 
appropriate choice of the time-dependence of i{}{t). Mod- 
els with log-normally distributed rates in discrete time 
have been considered in [J, [a] • 

Denoting the numeraire-rebased zero coupon bond 
prices as 
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FIG. 1. The solution of the model for the convexity- adjusted 
Libors Li for a simulation with n — 20 time steps (above), 
n — 40 (below), r = 0.25 and constant forward short rate 
ro — 5%, for several values of the volatility ip. 



defined by appropriate interpolation. However, for time- 
independent volatility "0, it can be shown that the model 
can be solved analytically. In this case the solution of the 
recursion equation for the one-step rebased zero coupon 
bonds Pi^i+i has the general form 



n—i—l 



P,^,+i{x,) = Y^ c} 



ii)p3->pXi-^{3i,fU 



(9) 



J=0 



„(») 



with Cj a set of constant coefficients. The convexity- 
adjusted Libors are given by 

f Po,i — ^0,1-1-1 



n—i — l 



(10) 






(«) 



The matrix of coefficients c^ has a triangular form 
(e.g. for n = 5) 
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(11) 



JS) 



The coefficients c, satisfy the recursion relation 



we note that the martingale condition for P^^j can be 
expressed as 



1 



^.,. =E[7^|^, 



(5) 



Two particular cases of this relation are 

P„^+i - E[A+i,,+2(l + Z,+ir,+i/,+i(a;))|J-,] (6) 

Po,. = E[Pm+i(1 + ^.^J»(2;))] (7) 

where we defined Ji(x) = exp(-!/;a;i — ^tp'^ti). 

These relations can be solved recursively for Pi^i+i and 
Li, starting with the initial conditions 



Ln~lTn-l — Po,n-l — 1 , Pn-l,n ~ 1 



(8) 



and proceeding backwards in time. In the next section 
we present a method for solving these recursion relations 
in analytical form. 



III. ANALYTICAL SOLUTION 

In practice the expectation values in Eqs. ([5]) are com- 
puted by numerical integration over the functional de- 
pendence of Pi.i+i on the Markovian driver Xi, which is 



„(») 



cf ^' + L, 



+ix,+i4!^)e(^-i)'^^*-+^ (12) 



which must be solved simultaneously with Eq. ()10p for Li. 
The initial condition is Cq — l,i„_iT„_i ~ Po,n~i — 
1. The recursion relation (IT^ can be solved backwards 
in time, for all i < n — 1, finding all coefficients in the 
matrix c starting from the upper left corner and going 
downwards. 

Once the coefficients c.^' and the convexity-adjusted 
Libors have been determined, all zero coupon bonds can 
be found as 



Pi,j[Xi) 



Pi,j(Xi) 



Pl,i+l{Xi)[l + LiTifi{Xi)] 



(13) 



where 



P^Ax^) = E[;^l-^.] = E[4.+i(l + L,r,e^^^-^^''^)\:F.] 



3,n 
n-j-1 



k ^ 



(14) 



fc=0 



n~j-l 



LjTj J2 



U)(k+l)if,x,-^{k^ + l)^^t,+k-4,^(tj-t,) 



cC e 



fc=0 



This completes the exact solution of the model. 



.(0 



The coefficients c^- satisfy certain general relations 
and sum rules, valid for arbitrary volatility. The first 
two coefficients Cq i can be given in closed form 



.(») 



Cq' — 1 for i = n — l,n — 2,- ■ ■ ,0 (15) 

n-l 



.(») 






J '3 



(16) 



A sum rule for the coefficients which will be useful in the 
following is 



E 

3=0 






(17) 



We illustrate the solution of the recursion in Figure [1] 
where we show results for the convexity adjusted Libors 
Li for several values of the Libor volatility tp- The two nu- 
merical examples considered assume a flat forward short 
rate tq ~ 5% and a time discretization with quarterly 
time steps r = 0.25. The terminal bond maturity is 
tn = 5 yr, and t„ = 10 yr, corresponding to n = 20,40 
steps, respectively. 

From Fig. [T] one observes that the convexity adjust- 
ment Lf"'^ — Li is always positive, and increases with the 
volatility ij). It is largest in the middle of the simulation 
interval, and with increasing volatility it becomes larger 
in a wider region expanding towards the beginning and 
the end of the simulation interval. This general behavior 
is expected on general grounds for the convexity adjusted 
rate Li provided that Li is positively correlated with the 



Libor associated with the payment delay Lt 



The 



convexity adjusted rate appears to vanish in the middle 
of the simulation interval, for sufficiently large volatil- 
ities. The rate of vanishing has a sudden increase for 
volatilities larger than a certain value. This phenomenon 
is the main subject of this paper, and will be explained 
and quantified in Section ITVl below. 



A. Scaling 

The model is uniquely defined by the parameters 
{U}, {Po,i}j V"- The zero coupon bonds Po,i can be equiv- 
alently written in terms of the zero rates r^ as Po.i = 
exp(-riii). 

The model is invariant under a simultaneous rescaling 
of these parameters given by 



ti — >■ Xti 

V'^ A-i/2^. 



(18) 



S^ 



Under this transformation, the coefficients c^ and the 
expectation values Ni are invariant, while the convexity 
adjusted Libors scale as Li — >• X~^Li. Due to this scahng 
invariance, the number of the relevant parameters of the 
model is reduced from three to two. 



B. Generating function 



In this section we present an efficient method for solv- 
ing the recursion relation p^ for the coefficients c^ . We 
introduce the generating function at the time horizon ti 

n — z — 1 



/«w- E 

J=0 



(») 3 

Cj X-' 



(19) 



This function satisfies a recursion relation, expressing the 
generating function at time ti in terms of the generating 
function at the next time tij^i 

/(') (x) = /('+!) (x) + Z,:+irx/('+i) {xe^"''^^ ) (20) 

The initial condition for the recursion is /'""^'(x) = 1. 
The expectation value Ni appearing in the expression for 
the convexity-adjusted Libor Li Eq. ^ is 



TV, = /W(eV't.) 



(21) 



The generating function f^^'{x) takes known values at 
x = 0, 1 



/W(0) = 1, /W(l) = Fo,»+ 



W( 



(22) 



where the second relation follows from the sum rule 

Eq. dni). 

In the zero volatility limit "0 = 0, the generating func- 
tion /(*) (x) can be found exactly 



fl,'\x)^n]-l^,{l + Lf\x). (23) 

where L^^'^ = t~^(Pq^j/ Pqj+i - 1) is the forward Libor 
rate for the time period (tj,tj+i). Expanding in powers 

of X this gives all the coefficients c, in the zero volatility 
limit. 

As the volatility increases ip > 0, the generating func- 
tion /'*'(x) also changes, in such a way that the con- 
straints Eqs. (P^ are still satisfied. 

The recursion relation for f^\x) can be reformulated 
in such a way that it does not contain any reference to 
the convexity-adjusted Libors 



/W(:r) = /(»+!) (:e) + (Fo,:+i - Po,^+2)x 



fii+i)^. 



^ru. 



/(«+i)(eV'"t.+i) 
(24) 



In the asymptotically large volatility limit ^/j — >■ oo, 
this recursion relation can be again solved exactly, and 
the generating function f^^' [x) is given by the asymptotic 
form 

/«(a;) = 1 + (Po.n-i - 1)2: + (Po,„-2 - Po.n-i)x^ 
+ ■ ■ ■ + {Po,.+i - Po,.+2)x"-'-' . (25) 

Note that the coefficients cf have well-defined limiting 
values as V' — >■ 00. 



This shows that the model considered has two very dif- 
ferent hmiting regimes, corresponding to i) small volatil- 
ity, and ii) large volatility. We will denote these regimes 
as the phases of the model. In each phase the generating 
function /(*' (x) has a well-defined expansion, given by 
Eqs. dSSl) and ^^, respectively. 

We can use the results for the generating function 
in the zero and large volatility limits in order to ob- 
tain asymptotic expressions for the convexity adjusted 
Libors Lj in the small and large volatility limits. This 
can be done using the relation between the expectation 
values Ni defined in Eq. pUj) . and the generating function 

We start by considering first the asymptotics of Li in 
the small volatility limit ij^'^ti <^ 1. The expansion of the 
generating function /q (x) around x = 1 reads 

ri-l rfwd 

/«(x) = Po.m[l+ E T7^^(^-l) (26) 

+0{{x-lf)]. 

Using Eq. (ITUl) this gives the small volatility asymptotics 
of the convexity-adjusted Libors L^, valid up to correc- 
tions of 0((^2^,)^) 

"-1 rfwd 

L.^L^l- E 1^71^7 (^"'*--l)) (27) 

ip^ti < 1 

(Note that the use of the zero volatility limit of the gener- 
ating function /q (x) was sufficient in order to derive this 
result. This is due to the exact condition /'•*^(1) = Po,i+i, 
which implies that the Taylor expansion of f^'^\x,i/j'^ti) 
in powers of a; — 1 and Tp'^U contains only a linear term 
in a; — 1 but not in ip'^U.) 

The low volatility approximation Eq. (|27)) has the 
familiar form of the convexity adjustment for a log- 
normally distributed rate. The growth of the convexity 
adjustment with the volatility has the familiar exponen- 
tial form, proportional to exp(V'^ti) — 1. Assuming a fiat 
forward Libor curve, we have L^^^^ = L^"'^ = i(e''o^ - 1). 
Then all terms in the sum over j are equal, and the result 
(P7| simplifies as 



U 



Lr(l-(n-z-l) 



L 



fwd- 



1 + P^<^Tj 



(e 



ip ti 



1))(28) 



This convexity adjustment is largest in the middle of the 
simulation interval, and vanishes near the boundaries. 
This agrees qualitatively with the main features of the 
convexity adjustment observed in Figure [T] 

Consider next the large volatility asymptotics i/'^ti ^ 1 
of the convexity adjusted Libors Li. This follows from 
the large x asymptotics of the large volatility generating 

function foo {x) 



Assuming again a flat forward Libor curve, this gives 
the large volatility asymptotics of the convexity adjusted 
Libors 



L, = L 



fwd 



1 + if"^T 



■'{n-i-l)ip^ti 



1p\ > 1 



(30) 



This shows that in the large volatility limit, the con- 
vexity adjusted rates Li drop off much faster with the 
volatility ip. The decrease is still exponential, but it is 
much faster due to the additional factor n — i — 1 (equal 
to the number of time steps to maturity) in the expo- 
nent Li ~ exp(— (ri — i — l)'0^ij). The vanishing of the 
convexity adjusted Libors is faster in the middle of the 
simulation interval, just like in the small volatility case. 
As an aside, we note that all dynamical quantities of 
the model can be expressed formally in terms of the gen- 
erating function f^^'{x). For example, the rebased bond 
prices are given by 



PM+i(a:)=exp(-lt,a2)/W(e'^-) 
P,,,{x)^eM-lud'.)f^'\e^n 



(31) 
(32) 



+LiTie'* 



4>x-^,P^U r(j) ^^lf,x+i,^(tj-ti)\ 



from which all zero coupon bond prices can be obtained 
using Eqs. ([Tc 



fg\x) ^ iPo,,+i - Po,,+2)x 



n~i—l 



0[x 



n—i — 2 \ 



.(29) 



IV. LIBOR VOLATILITY TRANSITION 

As the volatility increases from zero to a large value, 
the coefficients cj interpolate between the two limiting 
values, corresponding to low and asymptotically large 
volatilities, respectively. Equivalently, the generating 
function f^\x) changes between the two limiting ex- 
pressions Jq (x) and f6o (x), in such a way that the two 
constraints Eq. (P^ are still satisfied. 

As mentioned above, the convexity adjusted Libors Li 
appear to become vanishingly small at some value of the 
volatility, see Fig. [T] This phenomenon occurs first in the 
middle of the simulation time interval, and then it gradu- 
ally extends also towards the boundaries. This is related 
to the expectation values Ni which become very large as 
the volatility increases. To investigate this in more de- 
tail, we show in Fig. [5] the plots of log A'i as function of 
the volatility ip for two simulations. 

We observe that the change is not gradual, but hap- 
pens at a sharply defined value of the volatility, which 
will be called the critical volatility Tpcr- The transition 
becomes more sharp as n — i ~ 1 increases. The critical 
volatility demarcates two regions of very different quali- 
tative behaviour, in which the model has distinct limit- 
ing expressions for the functional dependence of discount 
bonds on the Markov driver. 
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FIG. 3. The zeros of the generating function f^''{x) in 
the complex plane for several values of the volatility ^ = 
0.5, 0.55, 0.6 around the critical value "i/icr = 0.53 (blue, vio- 
let, red). The black dots correspond to the large volatility 
asymptotic generating function fSc{x), which is obtained for 
ip = oo. The parameters are the same as in the upper plot of 
Fig. [2l and the time slice ti considered is i = 10 . The circle 
shown has radius exp{ip'^^ti) . 



FIG. 2. The log of Ni as a function of the volatility ip 
(solid line), and its large volatility approximation as Ni = 
f^\exp{Tp^ti)) (dashed line). Above: n = 20, i = 10. Below: 
n = 40, i = 30. Both cases correspond to r = 0.25, ro — 5%. 



The convexity-adjusted Libors Li become very small 
for volatilities above the critical value ip > V'cr- In prac- 
tice, they can become so small that they are below ma- 
chine precision. This phenomenon thus imposes a limit to 
the practical applicability of the model, and it is impor- 
tant to understand the conditions under which it occurs. 

In the following we investigate in some detail the na- 
ture of this transition, and formulate a criterion for find- 
ing the critical volatility V'cr at each time horizon ti . We 
will show that the singular behaviour of the expectation 
values Ni is related to the distribution of the zeros of the 
generating function f^'^'{z) in the complex plane. 

The generating function f''^\z) is a polynomial of z 
with positive coefficients, and thus does not have any 
zeros on the positive real axis z > 0. However, it is 
well known that the position of the zeros in the complex 
plane can influence the behaviour of the function along 
the real axis. The generating function /'■'■' (z) has n — 
i — 1 zeros. At zero volatility, all zeros are at large real 
negative values Zk — —l/[L^^^Tk), but they migrate in 
the complex plane as pairs of complex conjugate values 
as the volatility increases, and surround the origin, see 
Fig. [21 As the volatility increases to very large values, 
the zeros reach fixed positions, given by the zeros of the 



asymptotic generating function foo{x). 

Furthermore, as the polynomial order n — i — \ in- 
creases, the number of the zeros increases and they close 
in on the positive real axis, pinching it at some point 
z* > 1. We will show that the function /^^•'(z) is contin- 
uous at z^., but its derivative has a jump. Thus the gener- 
ating function has a cusp at this point. Recalling that Ni 
is related to the generating function as Ni = /(')(e'^ *•), 



see Eq. (|2ip . it follows that Ni has a singular behaviour 
at the volatility f/'cr, given by the equation 



exp(V'crii) = 2;* 



(33) 



This equation determines the critical volatility ^cr at the 
time slice ti. Geometrically, this has the following mean- 
ing: the critical volatility at the time horizon ti is given 
by that value of ip for which the zeros of the generating 
function f''^\z) enter the circle of radius exp(-0^ii). As 
ijj increases, the zeros move closer to the origin, while the 
circle of radius exp(?/'^ii) expands, such that at some in- 
termediate value ■0cr7 the zeros will cross the expanding 
circle. This picture is illustrated in Figure [3] on the ex- 
ample of the transition shown in the upper plot of Figure 

m 

The transition between the two volatility phases is sim- 
ilar to a phase transition in the Lee- Yang formalism [0] , 
where critical points are associated with the values of the 
thermodynamic parameter (fugacity) at which the zeros 
of the grand canonical partition function pinch the real 
axis. For a general introduction to phase transitions see 



Q. As shown in [9], the partition function is continuous 
at the critical point, but its derivative has a jump, which 
is proportional to the density of zeros around this point. 
These results hold in the thermodynamical limit of an 
infinite volume; the analog of the thermodynamical limit 
in our case isn — i— 1— >-oo, where n — i — l is the number 
of time steps from maturity to the time slice considered 
U. 

Such a behaviour is precisely what is observed in Fig.[2l 
where one can see that log Ni is continuous everywhere, 
but its derivative has a jump at ip„- The analog of the 
partition function in our case is the generating function 
f^^\z), and the jump occurs at z* = exp(i/)^jii) In the 
following we quantify this statement, and compute an 
explicit result for the discontinuity of the derivative of 
/W(z) at z,. 

The generating function can be written explicitly in 
terms of its zeros z^ as 



nrr^(i 



^fc=i 



z/zk) 



(34) 



As the polynomial order nf = n — i — 1 increases, the 
roots arrange themselves on a closed curve around the 
origin, which can be parameterized in polar coordinates 
as z{9) = p{9)e^^ . The roots appear in complex conjugate 
pairs, which implies that the curve describing the zeros is 
symmetric under reflection on the real axis p{9) — p{~9). 

The logarithm of the generating function f^'^''{z) can 
be expressed in the limit Uf -^ cx) as an integral 



log/«(z)=/ d9g{9)log{ 



pHO) 



2 cos 6* 



piO) 



K}35) 



where g{9) is the density of roots at polar angle 9. It is 
normalized as 



d09{9) = luf 



The derivative of log J*^*-' (z) is 

d 

dz 



log/«(z) = 2 / d9g{9) 



z — p{9) cos 9 



z^ - 2 p{9) cos 9z + p'^ 



(36) 



7) 



This is discontinuous across the boundary aX z — p(0) 
with a jump given by the density of zeros at this point 



dz 



log/W(^)| 



27rg(0) 



p(0)+e 
1 



dz 



logf'Kz)U=pio)- 



P(0) 



(38) 



These results are completely analogous to the expressions 
derived in f!l| for the jump of the derivative of the grand 
canonical function at a critical point. 

A similar phenomenon occurs for any expectation value 
of the form similar to Ni, with </> a real number 



m,^+ie 



(px-^cj?ti^ _ Hi)f^ip4>ti\ 



(39) 



The expectation value can be expressed in terms of the 
generating function f^^\x) as shown. The critical volatil- 
ity corresponding to this expectation value is found in 
analogy to Eq. (|33l) and is given by exp(i/'(/)ti) = z*. 
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FIG. 4. The integrand in the expression (|40p for the expec- 
tation value Ni at time ti — 2.5 for several values of the 
volatility ^p: i) 0.4 (solid line), 0.5 (dashed line), 0.52 (dotted 
line). The parameters of the simulation are n = 20, i — 10, 
ro = 5%, r = 0.25. Note the second maximum which appears 
for volatilities close to the critical volatility ^cr. 



We close this section with a comment about the prac- 
tical observability of the phase transition in usual imple- 
mentations of the Markov functional model. The expec- 
tation value Ni defined in ([TU| is written explicitly as 



iV,: 



E[A,.+i/,(x)] 
dxi 



(40) 



^/2^ 



'2*7 p. 



i,i+l 



i^XijC 



iixi-^ip^ti 



This integral is usually assumed to be dominated by 
contributions from small values of Xi around the origin 
la^il £ K^/ti, where k ~ 3 — 5. The integral is then evalu- 
ated numerically either on a grid, or using Monte Carlo 
simulations. However, for volatilities ip above the critical 
value, the integrand develops a second local maximum 
at large values of the Markov driver as seen in Figure SI 
which will dominate the integral above the critical volatil- 
ity. In the example of Figure 21 the secondary maximum 
appears at a; '-^ 12, which is almost 10 standard deviations 
away from zero. Thus the integral ()40p will receive signif- 
icant contributions from a region in x which is sampled 
very inefficiently in Monte Carlo or grid methods. This 
implies that the model will not be simulated correctly in 
the large volatility phase, and the phase transition will 
be unobservable under these simulation methods. 



V. EXAMPLE: CONSTANT FORWARD SHORT 
RATE 

We illustrate the general results discussed in the pre- 
vious section on the example of a forward yield curve 
with constant short rate rg. The initial yield curve is 
Po,i = exp(— roii), and the numeraire rebased discount 
bonds are Pqa = exp(ro(n — z)t). 



In Figure [3] we show the zeros z^ of the generating 
function f^'^\z) at the time shoe i — 10, for a simulation 
with ro = 5%, total simulation time tn = 5 yr, with time 
step T = 0.25 and n — 20 time steps. The colored dots 
show the exact zeros at three values of the volatility ip — 
0.5, 0.55, 0.6 around the critical volatility ipcr — 0.53, and 
the black dots are the zeros of the asymptotic generating 
function fH (z). 

We note that the asymptotic generating function 
/oo {z) gives a reasonably good approximation for the 
position of the zeros in the large volatility phase. For 
this reason we discuss in some detail the position of the 
zeros of the asymptotic generating function fsc (z) for 
which an analytical treatment is possible in the limit of 
a constant forward short rate. 

The asymptotic generating function for a constant for- 
ward short rate ro is 



/«(x) = l + [l 



a — i — 1 



'xY 



(41) 



Its zeros are Xk = e ''"^Zfe, where z^ are the zeros of the 
polynomial pn^ (z) of degree Uf = n — i — I 



Pnj {z) = — 



1 



+ z"/ 



(42) 



Under usual market conditions 1/(1 — exp(— rpr)) > 1, 
and by the Enestrom-Kakeya theorem [l3], all zeros of 
this polynomial lie outside the open unit disk \zk\ > 1- 

The structure of the zeros of p„j, {z) can be studied by 
noting that this polynomial is the truncated Taylor series 
of the function 



F{z) 



1 



1 



e''o^ - 1 1 - z 



(43) 



The function F{z) has an exact zero at zq = e''"'^, and a 
pole at z = 1, which means that the convergence region 
is the circle \z\ < 1. 

The theory of the zeros of a truncated Taylor expansion 
is a well studied subject in approximation theory The 
main result is the Jentzsch-Szego theorem [111 Il2l |. ac- 
cording to which the zeros of the truncated Taylor series 
of a function F{z) either converge to the zeros of i^(z) as 
n/ — > cxD, provided that they are inside its convergence 
region, or they accumulate on the boundary of the con- 
vergence region. For our case the latter situation applies, 
such that the zeros of p„ , (z) will accumulate uniformly 
on the \z\ = 1 circle as the polynomial order ny increases. 

This implies that the zeros of foo{x) close in on the 
positive real axis at x* = p = e"''"'^ as the polynomial 
order nj — n ~ i — 1 becomes large. This phenomenon 
is visible already at moderate values of Uf ^ 0(10), as 
seen from Fig. [31 

We can compute the properties of the generating func- 
tion at the critical point a;* by applying the general re- 
sults discussed above. Taking p{9) = p = e~^°'^ and 



g{9) = n//(27r), the integral in Eq. 
exactly with the result 



35|) can be performed 



log/^'(^ 



n/logp, 
Uf log X , 



X < p 
X > p 



(44) 



The jump of the derivative across the critical point 
X* = p is 



^log/«(. 



P + e)-—logfg\x 



at 
^5) 



In the large volatility phase the derivative of log f^ (x) 
is very large, and is of the same order of magnitude as 
expected from the asymptotic form Eq. ([^ . 

Finally, we consider the case of practical interest of 
fini te p olynomial order Uf = n — i ~ 1. It was observed 
in [IJ] that for the truncated Taylor series a^ -\- aiz + 
•••-(- a„z" of a function F{z), a good approximation for 
the moduli of the zeros \zk\ is obtained by neglecting all 
but the first and last terms (by considering the simpler 
polynomial p(z) = ao-l-a„z"), provided that the function 
F(z) does not have zeros within the convergence region. 
In our case of the polynomial Pn^iz), the corresponding 



polynomial p„ 



1/(1 - e"'^"^) -t- z^-f has zeros 



Zk = ( 



1 - e-''o 



T * 



(46) 



fc = l,2,--- 



1 



The critical volatility can be computed using equation 
Eq. ([33]) and is given by 



aroT + lll„ti 



1 



1 - e-''"^ 



xl/(n-j-l) 



(47) 



^i 



lin ■ 



l)r roT 



The minimum value of ipcY is reached at i — [n/2\ , in the 
middle of the simulation interval, where i(n — i — 1) is 
maximal. This agrees with the shape of the Li curves in 
Fig. [H where the critical volatility is first reached in the 
middle of the interval. Thus the practical applicability 
range of the model is restricted to volatilities smaller than 
the minimum critical volatility 



V'^ < (V'cr)min = 



1 



1 



n/2] 



■log( — )• 



roT 



(48) 



This expression satisfies the general scaling properties of 
the model Eqs. ([T8|) . 

In Table 1 we show the values of the maximally allowed 
volatility for several values of the total simulation tenor 
tn = nr, the time step r, and the short rate rg. The 
maximum allowed volatility decreases with the size of 
the time step r, with the tenor of the simulation i„, and 
with the short rate tq. 





f„ = 5 


fn = 10 


fn = 20 


tn = 30 


ro 


r = 0.25 


r = 0.5 


r = 0.25 


r = 0.5 


r = 0.25 


r = 0.5 


r = 0.25 


r = 0.5 


1% 


48.95% 


65.10% 


24.48% 


32.55% 


12.24% 


16.28% 


8.16% 


10.85% 


2% 


46.04% 


60.70% 


23.02% 


30.35% 


11.51% 


15.17% 


7.67% 


10.12% 


3% 


44.24% 


57.96% 


22.12% 


28.98% 


11.06% 


14.49% 


7.37% 


9.66% 


4% 


42.92% 


55.94% 


21.46% 


27.97% 


10.73% 


13.99% 


7.15% 


9.32% 


5% 


41.87% 


54.32% 


20.93% 


27.16% 


10.47% 


13.58% 


6.98% 


9.05% 



TABLE I. The maximal Libor volatility tp for which the model is everywhere below the critical volatility ipa, determined 
according to Ea. (l47p . for several choices of the total tenor i„, time step r and the level of the interest rates ro . 



VI. CONCLUSIONS 

We discussed in this paper the behaviour of a Markov 
functional model with discrete log-normally distributed 
Libors, as a function of the Libor volatility. The model 
can be solved exactly in the T— forward measure in the 
limit of a time-independent volatility, by a backwards 
recursion relation. Analytical results can be obtained for 
the functional dependence of all discount bonds on the 
Markov driver. 

We showed that the model has two volatility regimes, 
corresponding to small and large volatility, with very dif- 
ferent qualitative behaviour. As the Libor volatility in- 
creases, there is a transition between the two regimes, at 
an intermediate critical volatility. In the large volatility 
phase the convexity adjusted Libors are very small, and 
can be below machine precision. Thus the existence of 
the large volatility regime imposes a limit on the appli- 
cability of such a model, which is manifested as an upper 
bound on the allowed Libor volatility. 

We formulated the conditions under which this phase 
transition occurs, and showed that it is related to the po- 
sition of the zeros of an appropriately defined generating 
function in the complex plane. The transition appears 
for long simulation times, and small time discretization 
steps. A similar phenomenon is expected to occur also 
for the practically relevant but analytically more complex 
case of time-dependent volatility ip{ti). Also, the discus- 
sion of this paper is limited to the T— forward measure, 
and one expects that the details of the volatility depen- 
dence could change in a different measure, but not the ex- 
istence of a phase transition. We hope to report progress 
along these directions in future work. 



Appendix A: The zeros of the generating function 



We study here the distribution of the complex zeros of 
the asymptotic generating function /^^ (z) for a constant 
forward short rate tq. This is related to the problem of 
finding the zeros of the polynomial 



-0.5 




Pn{z) = a + z + Z^ + 



+ z'- 



(Al) 



-1.5 -1 -0.5 0.5 



FIG. 5. The zeros of the polynomial Pq{z) for a = 80. The 
solid line denotes the curve p{9) given in Eq. (|A2|) . and the 
dashed line shows the approximation for the moduli of the 
zeros |zfc| = a}'" obtained by retaining in Pn{z) only the first 
and last terms. 



where a is a real number larger than 1. They are the same 
as the zeros of the equation z"+-'^ + {a— l)z — a = from 
which z = 1 is excluded. Expressed in polar coordinates 
z — pe^^ one finds that the zeros are on the curve p{9) 
which is given by the solution to the equation 

^2(„+i) ^g,) ^ (^ _ -^^2^2 ^g,) _ 2a{a - l)p{e) cos 6 + 02) 

Figure [S] shows the exact roots of pg(z), along with the 
curve p{9) (solid line), and the circle with radius a^/" 
(dashed line), which was used in the main text. We note 
that the latter approximates the moduli of the zeros very 
well, as noted in |13| . 

The curve p{9) intersects the real axis at some point 
p(0). The solution of the equation (|A2[) at 9 = has 
always the solution p{0) = 1, and in addition it can have 
another 2 solutions, or none, depending on the values 
of n, a. The solutions different from 1 can be found as 
the solutions of the simpler equation h{p) = p"-~^^ — (a — 
l)p + a = 0. The function h{p) has a minimum at p* = 



[(a— l)/(n+l)]^/". We distinguish the 3 cases, according 
to the value of h{p^) 

l./i(p*) < ^ 2 solutions for p (A3) 

2.h{p*) = ^ p = p* (A4) 

3./i(/9*) > ^- no solutions for p (A5) 

These cases are obtained for n < [n^],n — [n^],n > [n^] 
respectively, where n* is the solution of the equation 



n n 



1 



n+l 



(A6) 



For a = 80 one has [n*] ~ 22. A good approximation for 
the solutions in case 1 is obtained by Taylor expanding 
h{p) around p* to quadratic order. This gives 



p* ± 



np\ 



n+l 



y \{n+ l)p" ^ ' 
In Fig. Owe show only the largest of the 3 solutions 



(A7) 



The angular distribution of the roots can be obtained 
from the imaginary part of the equation for z, which 
reads p" sin[(n + 1)9] = —(a — 1) sin0. For the subset of 
even n, this has real solutions for p only if is in one of 
the regions 



e ((2fc-l)- 



1' 



,2k- 



1 



), fc = 0, l,---n-^8) 



There is one solution in each of these intervals of equal 
angular opening, which means that in the large n limit, 
the angular distribution of the roots approaches a uni- 
form distribution g(9) = n/(27r). In the same limit, the 
curve p{0) approaches the unit circle p{0) = 1, as required 
by the Jentzsch-Szego theorem. 
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